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Abstract 

For the first time a three-dimensional (3D) chaotic billiard - the 3D Sinai 
billiard - was quantized, and high-precision spectra with thousands of eigen- 
values were calculated. We present here a semiclassical and statistical analysis 
of the spectra, and point out some of the features which are genuine conse- 
quences of the three dimensionality of this chaotic billiard. 
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Most of the work in the field of "quantum chaos" is related to systems with two degrees 
of freedom. This is a natural tendency, since two-dimensional (2D) systems provide the 
simplest, non-trivial examples of chaotic systems. In this respect, 2D billiards were subjected 
to thorough studies and served as good paradigms for the understanding of time-independent 
chaotic Hamiltonian systems. 

The main motivation to extend the study of quantum chaotic billiards to three dimen- 
sions, is to examine the semiclassical approximation: this approximation introduces an 
inherent error of order fi 2 to the energy levels (regardless of dimension) [1], while the mean 
energy spacing is proportional to h d } for systems with d degrees of freedom [2]. This might 
prohibit the resolution of single (high-lying) levels when d > 2. In the absence of a complete 
theory, it is important to test these ideas, to check whether the semiclassical approximation 
can be recovered by suitable /i-expansions [3,4], and to evaluate its usefulness for statistical 
purposes. Other motivation to study 3D systems is the existence of "real world" systems 
that cannot be reduced into 2D idealizations, and which can exhibit some (classical and/or 
quantum) genuine 3D effects. Also, for d > 2, classical systems can support all the possible 
types of stability, including mixed (lexodromic) stability, which is absent in 2D systems. 

In this Letter, we would like to report about the numerical quantization and spectral 
analysis of the 3D Sinai billiard. The results will be presented with special attention to gen- 
uine 3D features which affect the spectrum in whatever way it is approached and analyzed. 
To the best of our knowledge, this is the first study of 3D quantum chaotic billiards which 
can match the studies of 2D billiards in accuracy and extent. We study 3D Sinai billiard, 
which is the space inside a cube of edge L and outside an embedded sphere of radius i?, 
R < L/2 (see fig. 1). The classical motion consists of a particle that is specularly reflected 
from the sphere and from the the cube's walls. 

The 3D Sinai billiard was quantized using the KKR method [5-7] which was previously 
used for the quantization of the 2D Sinai billiard [8]. The essence of the method is to 
transform the Lippmann-Schwinger integral equation into a particularly efficient system of 
linear equations. This is obtained by a special choice of the basis functions and exploiting 
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the symmetries of the problem: the lattice periodicity and the spherical symmetry of the 
inscribed potential. This reduces the numerical effort considerably. The method heavily 
relies on a precise and efficient calculation of "lattice constants" [6] which is achieved by 
using the full Ewald summation technique [7,9], and choosing the convergence parameter in 
an optimal way which leads to the fastest convergence for a prescribed accuracy. Details of 
the calculation will be given elsewhere [10]. 

The Hamiltonian of the full 3D Sinai billiard is invariant under the discrete group Oh of 
the cube's rotations and inversion [11]. This group consists of 48 elements and 10 classes, thus 
the quantum system suffers from geometrical degeneracy, and the spectrum is composed of 10 
independent spectra. In order to remove this "trivial" degeneracy [8], we desymmetrized the 
billiard by using as basis functions linear combinations of spherical-harmonics that transform 
according to one of the irreducible representations (irreps) of Oh ("cubic harmonics" [12]). 
We found a simple, non-recursive method for determining the cubic harmonics, using the 
singular value decomposition (SVD) algorithm [13], and calculated the coefficients for the 
completely symmetric (Al) and completely antisymmetric (AY) one-dimensional irreps. 
These correspond to introducing Neumann and Dirichlet boundary conditions, respectively, 
on the symmetry planes of the cube, forming a desymmetrized element which is a triangular 
pyramid (it is marked by heavy lines in fig. 1). It is important to emphasize, that the above 
decomposition of the Hamiltonian is crucial for the feasibility of the calculation: had the Oh 
symmetry not exist, we would have to deal with matrices larger by a factor of 48 2 . 

We calculated few independent spectra of the billiard for L = 1 and R = 0.2, R = 0.3. 
Dirichlet and Neumann boundary conditions were applied independently on the sphere and 
on the symmetry planes. The spectrum for R = 0.2 and Dirichlet boundary conditions on 
both sphere and planes contained 2288 levels and covered the range < k n < 200 (k 2 n = E n ). 
This is the largest spectral interval obtained so far, and we shall refer to it as R = 0.2 (DD). 
The other three spectra for the same radius contained 1000 levels each, and the spectra for 
R = 0.3 contained 760 levels each. The completeness of the data was checked by comparing 
the actual number of eigenvalues N(k) to Weyl's law N(k) [14], in order to detect either 
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missing or redundant eigenvalues [15]. For all the spectra computed, N(k) — N(k) fluctuated 
about with no systematic deviations, indicating with high probability the completeness 
of the spectra. A further consistency check of the spectrum is obtained by comparing the 
lower eigenvalues to the eigenvalues of the R = integrable (and separable) billiard. For 
Dirichlet boundary conditions on the planes, the lowest angular momentum that appears in 
the KKR determinant is / = 9. Thus, for kR < 9 we expect the eigenvalues of the billiard to 
be very close to the R = eigenvalues, which is indeed confirmed by the numerical results. 
We performed also some convergence checks to estimate the error in the eigenvalues, and 
found it to be of order 10~ 4 of the mean level spacing. 

In order to facilitate later analysis, it is instructive to examine first the various classes 
of the classical periodic manifolds that exist in the 3D Sinai billiard, since they constitute 
the building blocks of the semiclassical analysis. First, we have the generic, isolated and 
unstable periodic orbits (1-dimensional manifolds). They are typical to hyperbolic systems 
in any dimension, and their total number proliferates exponentially with their length [2]. For 
our billiard, they are characterized by at least one non-tangential collision with the sphere. 
Second, we have the non-generic, neutrally stable, "bouncing ball" orbits [16,8,17,18], which 
bounce only between symmetry planes and/or the walls of the cube. For the specific case 
of the desymmetrized 3D Sinai billiard, they are conveniently analyzed by considering the 
integrable R = case. The quantum density of states can be explicitly written in this case 
as: 

d R=0 (E) = J2 8 

n>l>m>0 




where sinc(x) = sin(x)/x, and p } q } r £ Z. The oscillatory part of the first term 
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E ~ (f )V + /2+m2 ) 



((p } q } r) ^ (0,0,0)) describes contributions from periodic manifolds (tori) that occupy 3D 
volumes in configuration space. Each such manifold is a (2-parameter) continuous family 
of periodic orbits of the same length. The other oscillatory contributions come from fD 
and 2D manifolds of periodic orbits that are confined to special lines and planes, respec- 
tively. They involve collisions with corners and edges, and need special care which will be 
given elsewhere [10]. It is interesting to note, that isolated and neutral periodic orbits (ID 
manifolds) exist, and in fig. 1 we show two of them, that contribute to the fifth and sixth 
terms of (1). The triangular face ABC supports an infinite set of 2D manifolds that stay 
on it, as if it were a 2D triangle. The number of bouncing ball manifolds grows as a power 
of their length, in sharp contrast to the exponential growth of the generic unstable orbits. 
For R > 0, infinitely many periodic manifolds are either completely or partially pruned. 
The remaining periodic manifolds are "terminated" with periodic orbits that are tangent to 
the sphere. These tangential manifolds constitute the third kind of periodic manifolds, and 
they are of one dimension lower than the bouncing ball manifold that they terminate. A 
comparison with the situation in 2D Sinai billiard reveals a substantial difference: there is 
only a finite number of bouncing ball manifolds in the 2D case (0 < R < L/2) [8], in contrast 
to the infinite number in the 3D case. Also, only in 3D case there exist 3D bouncing ball 
manifolds and 2D tangential manifolds. These differences emphasize the almost unavoidable 
role of the bouncing balls in the 3D case, which is a direct consequence of the high symmetry 
of the original, non-desymmetrized billiard. Thus, the same high symmetry that allowed 
an efficient decomposition of the Hamiltonian, is responsible for the prominent non-generic 
bouncing ball manifolds. 

In the sequel, we shall study the relations between the quantum spectra and the classical 
periodic manifolds, in light of the discussion above. The first step is to construct the length 
spectrum, which is the (finite) cosine-transform of the oscillatory (wavenumber) density of 
states: 
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Semiclassically, D(x; k) is expected to be peaked near lengths x that correspond to classical 
periodic manifolds of the billiard [16,17], thus, allowing us to extract classical information 
out of the quantum data. Substituting (1) into (2), we get that the contribution of a d- 
dimensional manifold to D(x; k) is 0(k^ d ~ 1 ^ 2 ). This is also valid for the ID unstable periodic 
orbits [16]. Hence, the (leading, 3D) bouncing ball contributions are expected to be much 
stronger (0(k)) than those of the isolated periodic orbits. This is to be contrasted with 
the 0{k 1 / 2 ) enhancement of bouncing balls in the 2D case [16]. The length spectrum of the 
R = 0.2 (DD) spectrum is plotted in fig. 2 for k = 100 and k = 200, and the lengths of the 
bouncing ball orbits for < x < 2 are marked together with their expected ^-dependence. 
We clearly observe peaks for the predicted lengths (with possible exception of ID contribu- 
tion at x = 2/y/3 ~ 1.15 due to interferences), and qualitative agreement with the predicted 
^-dependence. The bouncing balls are clearly seen to dominate D(x; k). Our only hope to 
detect an unstable periodic orbit is to look in the domain where the contributions are still 
fairly isolated, x ^ 0.6. The shortest unstable periodic orbit, along the edge AE of fig. 1, 
has length of 0.6 for the R = 0.2 case, which is very close to the length l/y/3 ~ 0.577 of the 
shortest ID bouncing ball. For R = 0.3 the situation is better, because the shortest unstable 
periodic orbit has length of 0.4, which is well separated from the shortest ID bouncing ball. 
In both cases, however, the expected contribution is not observed. A detailed calculation of 
the semiclassical contribution of this unstable periodic orbit reveals, that being on the axis 
of an 8-fold symmetry greatly diminishes the semiclassical contribution of this orbit with 
respect to the situation where desymmetrization is not imposed. Quantitatively, the ratio 
is given by: 



M = j 



2a/1 -2/3-/3 



/3 N 



1-/3, 



« (^j for -> 

In particular, p 8 (R/L = 0.2) 2 X 10~ 4 <C 1. We emphasize that this is a purely 3D 
geometrical effect due to desymmetrization, and that the amplitude of this contribution to 
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D(x; k) is ^-independent as for all other unstable periodic orbits. 

The above analysis shows, that if one is interested in the study of the generic contribu- 
tions, it is imperative to develop a method which eliminates the non-generic contributions 
in an accurate and simple way. One possible idea is to parameterize the boundary condi- 
tions on the sphere, and use this freedom as a selection filter to extract the contributions of 
generic periodic orbits exclusively. Work along these lines is now in progress [19]. 

We applied some of the common spectral statistics to the spectra described above. In 
fig. 3 we present the nearest neighbor distribution P(s) for the R = 0.2 and R = 0.3 spectra. 
Both curves are similar to the GOE curve, but for R = 0.2 there are noticeable deviations 
for small s. These deviations are attributed to the bouncing ball orbits described above. 
The better agreement for R = 0.3 is qualitatively explained by the pruning of bouncing ball 
manifolds and the reduction of their measure as R increases. 

As a representative of the bilinear statistics (which are based on the two-point correlation 
of the spectral density) we show E 2 (7) [20] for R = 0, 0.2, 0.3 and Dirichlet boundary 
conditions on the planes (fig. 4). The energy interval for which S 2 (/) is calculated is taken 
to be 133 < k < 146, such that d(E) does not change significantly within this interval. 
For the chaotic cases, R = 0.2, 0.3, the departures from the universal GOE curve occur 
earlier than predicted by Berry [21]. Also, in contrast with generic chaotic systems, S 2 (/) 
oscillates around a saturation value which is larger than the GOE value. These deviations 
are more pronounced for R = 0.2 than for R = 0.3, which is a clear indication of the effect 
of the non-generic bouncing ball manifolds, that are larger both in number and in measure 
in the R = 0.2 case. The bouncing ball manifolds attain their full strength in the R = 
(integrable) case, for which we observe an initial slope which is much steeper than the value 
1 predicted for Poisson statistics. This is due to number-theoretic (NT) degeneracy of the 
energy levels, which depends on the dimensionality. Taking this degeneracy into account, 
the initial slope of S 2 (/) for R = was calculated to be 0.0323 kL } which agrees well with the 
quantum data. (This estimate is based upon a numerical evaluation of the degeneracy, and 
was checked up to kL = (9(10 3 ).) The expression of the leading term in the ^-expansion of 
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the saturation value (R = 0) is semiclassically estimated by 1.36 10~ 4 (kL) 2 (\n(kL) — 0.96), 
which is due to both the large contributions of the bouncing balls in three dimensions, and 
a number theoretical degeneracy of their lengths. This saturation value was found to agree 
with the quantum results for very large k. For the case shown in fig. 4, k is not large enough, 
and the asymptotic value is not reached. 

To summarize, the quantization of the 3D Sinai billiard using the KKR method became 
feasible due to the high symmetry of the billiard, and the resulting block-diagonal structure 
of the Hamiltonian. This high symmetry results in non-generic bouncing ball periodic mani- 
folds. The infinite number of the bouncing ball manifolds were shown to dominate the length 
spectrum in a way that have important qualitative differences from the 2D case. Genuine 
3D geometrical effects give rise to dramatic reduction of the semiclassical contribution of 
the shortest unstable periodic orbit. The quantum short range level statistics were shown to 
agree quite well with GOE results, in accordance with the common experience with systems 
that have time-reversal symmetry. However, important deviations have been observed in 
the E 2 statistics which were attributed to the bouncing ball manifolds. 
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FIGURES 
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FIG. 1. Desymmetrized 3D Sinai billiard. Fhe billiard is indicated by boldface edges, and 
its corners are marked by letters. Chained line: shortest isolated periodic orbit of length Z/a/3. 
Double chained: isolated periodic orbit of length L/\/2. 
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FIG. 2. Absolute value of length spectrum for the R = 0.2 (DD) spectrum. Full line: k = 100, 
dashed line: k = 200. The first 8 bouncing ball orbits are indicated together with their relative 
contribution to D(x;k). 
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FIG. 3. Nearest neighbor distribution P(s). For each R, data are averaged over available 
spectra with different boundary conditions. Inset: difference between P(s) and Pgoe( s )- 
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